Giant anharmonicity suppresses superconductivity in AIH3 under pressure 
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The anharmonic self energy of two zone boundary phonons were computed to lowest order for 
AIH3 in the PmSn structure at 110 GPa. The wavevector and branch index corresponding to these 
modes are situated in a region of phase space providing most of the electron-phonon coupling. The 
self energies are found to be very large and the anharmonic contribution to the linewidth of one of 
the modes studied could be distinguished from the electron-phonon linewidth. It is found that an- 
harmonicity suppresses the electron-phonon coupling parameter A, providing a possible explanation 
for the disagreement between experiment and previous theoretical studies of superconductivity in 
this system. 



I. INTRODUCTION 

It has been suggested four decades ago that elemen- 
tal hydrogen could form an exceptionally high Tc super- 
conductor under compression^ (a recent estimate being 
242 K at 450 GPai^), and more recently that it could have 
very exotic properties, such as being a metallic quan- 
tum liquid^, or forming protonic Cooper pairs^^^. How- 
ever, metallic hydrogen has long been elusive; the dimers 
persist^ and hydrogen remains non-metallic up to a pres- 
sure of 320 GPa^ 

An alternative route to hydrogen superconductivity 
has been suggested in the form of hydrides, where the 
presence of a heavier element can act to chemically "pre- 
compress" hydrogen, compelling it to reveal its super- 
conducting properties at lower pressures^i^. Because 
of their large hydrogen content, superconductivity in 
the group IV hydrides has been studied extensively^Tii^. 
Some tri-hydrides also received attention recently^^ and 
it was suggested that the origin of superconductivity 
in these systems could be soft phonons in the vicinity 
of phase transitions. Within this context, aluminum 
hydride (AIH3) under pressure has recently been stud- 
ied both theoreticallysii — and experimentally^^. Us- 
ing random structure searching, a particularly interesting 
phase of symmetry Pm3n has been found to be energet- 
ically favorable above ^ 70 GFdJ^^. This phase con- 
tains two formula units per cell, with the Al ions forming 
a body-centered-cubic (bcc) structure and the hydrogen 
ions forming linear chains on the faces of the cubic cell. 
Ah initio calculations also suggested that the electron- 
phonon coupling parameter should be fairly large in this 
phase at ~ 110 GPa (A 0.74), predicting a value of 
Tc— 24 K,^^ in agreement with the general idea that com- 
pressed hydrides could be good superconductors'". Very 
interestingly, however, no superconducting transition was 
found down to 4K^^; the reason for the disagreement be- 
tween theory and experiment is unclear. 

In the present work, we show that the phonon modes 



which provide most of the electron-phonon coupling are 
actually strongly renormalized by anharmonicity, which 
should greatly affect the value of the predicted Tc. 

In sections |TT1 and IIIH basic formulas pertaining to 
phonon mediated superconductivity and phonon anhar- 
monicity are reminded, which also serves to fix the nota- 
tion. The ah initio calculations performed are described 
in section ITV} and the main results pertaining to anhar- 
monicity are presented in section FVl 



II. SUPERCONDUCTIVITY 

The theory of phonon-mediated superconductivity is 
well understood^^. A popular approximation to the su- 
perconducting transition temperature is given by the 
Allen-Dynes modification of the McMillan formula, ^^^^^ 
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(1) 



where /i* is a parameter of order 0.1 which approximately 
accounts for the electron-electron repulsion at the Fermi 
level (which tends to weaken Cooper pairs, and thus re- 
duce Tc), oJiog is the logarithmic average of the phonon 
frequencies and A is the electron-phonon interaction (or 
electronic mass enhancement) parameter. This last pa- 
rameter in turn is obtained from the Eliashberg spectral 
function a^F, 



A = / dujX{uj)] X{uj) = 
Jo 



a^F{uj) 



(2) 



The EHashberg spectral function can be approximately 
related to the phonon linewidths 7iy(q) by^ 



a^F{uj) 



1 



2iThN{eF) N ^ cj^(q) 



5{uj -ujj,{q)), (3) 



where TV is the number of unit cells in the crystal, N^cf) 
is the density of states per unit cell at the Fermi energy. 
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q is a wave vector constrained to the first Brillouin zone 
(IBZ), u is di mode label and cjj^(q) is the frequency of 
phonon mode qu. This implies that A can also be ex- 
pressed as 



with 



q,i/ 



with 



(4) 



(5) 



The usual method employed to obtain A from ab initio 
calculations is to first obtain the band structure of the 
system, second to obtain the phonon frequencies within 
the Born-Oppenheimer approximation, and third to ob- 
tain effective electron-phonon coupling parameters. The 
system is then approximately described in terms of a 
Frohlich Hamiltonian, 

nkcr c^u 



qi^ m,n\<.(T 



^qi^ )^m1k+q,nk' 



the parameters of which are set to the ab initio computed 
values. In the above m, n are band labels, a is a spin la- 
bel, k, q are wave vectors in the IBZ, {c^c'^} and {b/b^} 
are electron and phonon ladder operators, {e} are elec- 
tronic eigenvalues and {g} describe the strength of the 
scattering between electrons and phonons. It is standard 
to set the electronic energies to the Kohn-Sham eigenval- 
ues, the phonon frequencies to the Born-Oppenheimer 
frequencies and to extract the values of the g parame- 
ters from the deformation potential. Standard field the- 
ory methods^^ are then employed to derive the phonon 
linewidths, 7^^, from this Hamiltonian. 



III. PHONON ANHARMONICITY 

The position operator for the ions in a crystal can be 
represented as 



r,(R) = R + b, + u^(R), 



(7) 



where R is a lattice vector, is the basis vector for ion 
and u,^(R) is the operator representing the displacement 
of the ion from its equilibrium position. Within the adi- 
abatic approximation, which assumes the electronic sys- 
tem instantaneously adapts to the ionic positions, the 
total energy as a function of the ionic positions can be 
taken as an effective potential for the ions which thus 
dictates their dynamics. This potential is expressed as 



U[{u}]=Uo + J2^n[{u}], 



(8) 



x$^J:--"(Ri,...,R„), (9) 

where the Greek symbols ai,...,^^ represent cartesian 
coordinates. It is assumed that the crystal is stable and 
that, consequently, the linear term in the displacements 
vanishes identically. The dynamics of the ionic degrees of 
freedom are then described by the effective Hamiltonian 



(10) 



where T is the kinetic energy operator of the ions. It is 
convenient to consider a canonical change of variable to 
reciprocal space of the form 



(11) 



In terms of these new position-like variables, the poten- 
tial terms can be expressed as 



Un[{u}] 



n\ N 



{a/tq} 

x*"j:::":(-qi'-'-qJ; (12) 

some useful symmetry relations pertaining to these an- 
harmonic coefficients are reminded in appendix [Al 



A. Harmonic phonons 

If the potential energy expansion is truncated after the 
second order, the resulting approximate Hamiltonian is 
harmonic, and leads to the standard small oscillations 
problem. In this case: 

^2 = ^ E E [<;(q)]^^:::S(q)<l(q)' (13) 

{an} q 

where the usual dynamical matrix has been defined as 

z)X'(q)^$^j^,nq,-q)- (14) 

It is standard to consider a canonical transformation to 
ladder operators of the form 

K (q) = Y ^Zv (q)^qi^ ; = 



bl^,. (15) 



The displacement vectors are defined as 

:E«<.(q) 



^(q) 



h 



2M«w^(q) 



(16) 



n=2 



where is the mass of ion ti and the polarization vec- 
tors E, which are chosen to be orthonormal, are solutions 
of the hermitian eigenvalue problem 

^ /^fi^ E-iM) = u^MfE^iMl (17) 
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which also yields the harmonic phonon frequencies. It is 
useful to define a "mode mass" M., as 



M,(q) = 



EJx..(q)P 



(18) 



this quantity is then a gauge of what type of ions are 
involved in a given mode. The harmonic part of the ionic 
Hamiltonian can finally be expressed as 



where 



and 



(19) 
(20) 



B. Phonon-phonon interaction 



When the anharmonic contributions to the poten- 
tial cannot be neglected, the ionic Hamiltonian can be 
expressed in terms of the harmonic Hamiltonian plus 
phonon-phonon interaction terms. The anharmonic coef- 
ficients can be expressed in terms of the harmonic basis, 
and lowest order contributions to the self-energy anhar- 
monic correction for the mode (qy) are given by^ 



i^i,qi \ / 

U\f^{ci,uj) = -— 2^2^$,,,,,,,,(-qi,qi,0)^,,,,,,(0,q,-q) ' , 

^i^H^^^) = Yl XlXl^qi+q2+q,G|^^,^i,^2(q,qi,q2)P^(^,^^i(qi),^^2(q2)), 



(21) 

(22) 
(23) 



qi,q2Z^iz^2 G 



+ UJ2^ (^1 + nsiuji) + nB(c^2)) ^(^^i - ^2) (nB{uj2) - nB{uJi)^ 



^^,...^jqi,...,qn) = Y ^^i':::^:(qi'-'qnX^i(-qi)-<:^J-qn). 

{aK} 



(24) 



(25) 



In the above, the quantity ub refers to the usual bosonic 
occupation factor. The labels (T), (L) and {B) refer 
to the "tadpole", "loop" and "bubble" diagrams, as 
schematically represented in Figure [H It is noteworthy 
that only the "bubble" contribution actually depends on 
the frequency cj, and that only this term will have an 
imaginary contribution. Furthermore, the "tadpole" di- 
agram vanishes by symmetry in this system. The results 
above will be specialized to the point qx = 7r/a(0,0, 1) 
on the side of the zone, at zero temperature. 



IV. COMPUTATIONAL DETAILS 

Electronic properties were computed using den- 
sity functional theory (DFT) as implemented in the 
QuANTUM-ESPRESSO packagers. The exchange- 
correlation was treated using the generalized gradient 
approximation (GGA) of Perdew, Burke and Ernzer- 
hof (PBE)2ii2^. Ultrasoft pseudopotential^ were used, 
where 3s and 3p states of aluminium were treated as 
valence. The plane- wave basis cutoff was set to 80 Ry. 




(a) Loop (b)Bubble (c)Tadpole 



FIG. 1. Lowest order anharmonic self-energy diagrams. The 
lines represent phonon propagators and the vertices third or- 
der ("bubble" and "tadpole") and fourth order ("loop") an- 
harmonic coupling. 



First Brillouin Zone (IBZ) integrations were performed 
as sums on a 24 x 24 x 24 Monkhorst-Pack k-mesh, us- 
ing a smearing parameter of 20mRy. Phonon proper- 
ties were computed using density-functional perturbation 
theory (DFPT)^i^. Interatomic force constants (IFC) 
were obtained from dynamical matrices computed on a 
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12x12x12 q-mesh. The electron-phonon coupling com- 
putations were performed using electronic and phonon 
quantities interpolated on a fine 72 x 72 x 72 mesh. 



V. RESULTS 

The phonon dispersion was calculated for the PmSn 
structure with a lattice constant a = 5.82 ao, which 
yielded a computed pressure of 109 GPa. The phonon 
spectral function a^F, as well as the electron-phonon 
coupling parameter A were also computed; results can 
be seen in Fig. [2l The values we have obtained are 
A 0.61 and hujiog ^ 68meV, leading to 12K < T^ 
< 19K (0.14 > /i* > 0.1). As is clear from Fig. H 
the bulk of the contribution to A comes from narrow re- 
gions in the IBZ centered at X, for modes at ~ 20meV 
and ~ 85meV; attention has been focused on the rele- 
vant modes at X, assuming that they are representative 
of modes in that region of the IBZ. These doubly de- 
generate modes will henceforth be referred to as Xi (20 
meV) and X2 (85 meV) and will be labelled as ux- In- 
terestingly, for these modes the motion of hydrogen ions 
is perpendicular to their chains; the displacements for X2 
are mostly that of hydrogen ions (with a mode mass of 
1.1 Mh)^ whereas the displacements of Xi involve both 
types of ions (mode mass of 5.7 Mh)- 

To better understand the origin of the large contribu- 
tion to A of modes around X, frozen phonon-perturbed 
bands were computed for Xi and X2. A displacement of 
the form 

ujR;7?) = 7/e^^^-^x,,^(qx) (26) 

was imposed on the ions in a supercell geometry for var- 
ious values of the unitless parameter^^ r] (the displace- 
ment can be chosen to be real because of the symmetry 
of the qx point). The appropriate supercell corresponds 
to a doubling of the original cell along the c axis, yielding 
a tetragonal cell. In the presence of this frozen-in pertur- 
bation, the Kohn Sham Hamiltonian can be expressed to 
first order in r] as 

flKS ^ Yl [^nkcj^ka^nka (27) 
nkcr 

m 

The correction to the eigenenergies e^k will be of order 
77^ at a generic k point. However, in the case of band 
degeneracy, the usual response formalism breaks down 
and the change in the band energy can be linear in r]. 
The unperturbed cubic system has Fermi sheets centered 
at the R (doubly degenerate) and M (non degenerate) 
points. As can be seen in Fig. [3l these two Fermi sheets 
are centered about the M point in the supercell geometry. 
The two sets of bands intersect close to the Fermi energy 
along the F-M direction at a point ko; at this point, the 



Kohn Sham Hamiltonian can be modeled, to linear order 
in 7^, as 

/ eo r]g\ 
hKS= \ eo V9 ] . (28) 

where 

eo = Cnko; 9 = ^nio+qx,nko- (^^) 

The band label indicates one of the three r] = degen- 
erate states, and it is assumed that the only relevant 
coupling to linear order is between the non-degenerate 
M band and the doubly degenerate R bands. The eigen- 
values of this Hamiltonian matrix are given by 

e = 0,±V2rj\g\ (30) 

and it is straightforward to extract a value for \g\. We 
find that \g\ 440 meV for Xi and \g\ 470 meV for 
X2 . These couplings are very large on the phonon energy 
scale, and suggest that the origin of the large linewidths 
lies in strong scattering between the Fermi sheets men- 
tioned above. The large coupling associated to Xi and 
X2 prompted further investigation of these modes. 

In order to gauge the anharmonicity of these modes, 
total energy frozen phonon calculations were also per- 
formed. A displacement of the form given by Eq. (|26|) 
corresponds to 

<(q;r,) = r,v^<,^ (qx) W- (31) 

It is important to note that, for a finite value of 77, this 
does not correspond to a realistic configuration of the 
ions. Indeed, the zero-point energy of a mode is in the 
order of meV, an energy that must be shared by all ions 
(a number of order N). Thus, for a single mode, only 
an infinitesimal amount of energy can be assigned to any 
ion, leading to infinitesimal average displacement. It is 
then the sum on all modes that yield a finite average 
displacement for any ion. This discussion does not in- 
validate the frozen phonon calculations, as they are only 
performed to extract anharmonic coefficients. 

According to the expressions for the anharmonic en- 
ergy as a function of displacement, Eq. (j8j) and ([12]), the 
total energy per unit cell for a given displacement should 
then be 

^ = ^ + ^fi-.. (qx) + |J*4... + Oirj') (32) 
where 

^X (qx,qx,qx,qx) (33) 

has been defined for convenience (note that qx and — q^ 
are equivalent points by reciprocal lattice periodicity). It 
is straightforward to extract the values of ^cc;jy^(qx) and 
^A,vx from the energy as a function of r] using a finite 
difference scheme; results can be seen in figure |4j As 
can be seen on the figure, the quartic contribution to the 
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FIG. 2. (Color online) (Left panel) Phonon dispersion relation computed using DFPT. The areas of the red circles overlapping 
the phonon dispersion are proportional to Az^q. It is clear that the bulk of the electron phonon coupling can be attributed to 
modes in a region near X with frequencies ^20meV(Xi) and 85 meV(X2). (Right panel) Eliashberg spectral function and 
phonon density of states PDOS (both in arbitrary units), as well as partially integrated value of A. 



potential is very large, with <l>4 c:^ 155 meV for Xi and 
$4 223 meV for X2. In the naive approximation that 
the quartic anharmonic coupling is constant throughout 
the zone and that coupling to other modes can be ne- 
glected, namely 

*^»^x,^x,i^,i^(qx,qx,q, -q) ^ 

(qx,qx,qx,qx), (34) 

and at zero temperature, this leads to a frequency renor- 
malization through the " loop" diagram of 77 meV for Xi 
and lllmeV for X2 (see equation [2T]) . 

Drawing conclusions from results at a single q point 
can be premature, as is examplified by the case of MgB2: 
frozen phonon calculations similar to those presented 
above suggested that the E2g modes at F should be 
highly anharmonic^^; refined calculations of the "loop" 
and "bubble" diagrams for this mode revealed that in 
fact this is not the case^i^. The point is that a frozen 
phonon calculation provides no information on the cubic 
coupling, which can largely cancel the quartic contribu- 
tion, nor does it account for the fact that the anharmonic 
coefficients have dispersions {ie are functions of q). Nev- 
ertheless, the unusually large value of ^4 prompted a 
more thorough investigation of anharmonicity for these 





fUjOjy{Cix) 


^^,-x 




FP DFPT 


FP SFD 




19.5 19.4 


154.6 154.9 


X2 


86.9 86.8 


222.8 224.2 



TABLE I. Comparing harmonic frequencies and $4 parame- 
ters (in meV) obtained by different methods. DFPT means 
"density functional perturbation theory", FP means "frozen 
phonon" and SFD means "supercell finite difference" . Results 
are seen to agree very well between different calculations, giv- 
ing confidence on the convergence of the SFD method. 



modes. 

The necessary coefficients were obtained by finite dif- 
ferencing of dynamical matrices computed for appro- 
priate supercells (henceforth the "supercell finite differ- 
ence", of SFD, method); a complete discussion of the 
formalism can be found in appendix [B1 Briefly, anhar- 
monic parameters were obtained by computing dynami- 
cal matrices with ions slightly displaced according to the 
polarizations of modes Xi and X2. A centered, five points 
finite difference scheme was applied to these dynamical 
matrices, yielding partially mode-projected anharmonic 
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FIG. 3. (Color online) Electronic bands near the Fermi energy 
in the tetragonal supercell. The bands for the unperturbed 
structure are shown (full black line), as well as for a X2 frozen 
phonon displacement corresponding to 77 0.4 (dashed red 
line). The zero of energy is set at the Fermi energy of the 
unperturbed system. The Xi data is not significantly differ- 
ent. The dashed vertical line indicates the point ko where 
the unperturbed bands cross. (Inset) Perturbed eigenvalues 
as a function of 77 corresponding to the degenerate energy in- 
dicated by the dashed line. The splitting is linear with the 
distortion for two bands (circles and triangles) and remains 
quadratic for the third (squares); the value of |^| : is directly 
related to the slope of the change in band energy with 77, and 
is equal to ^ 470 meV in this case. 









n. 






shift (%) 




22.9 


-6.3 


16.6 


19.4 


31.9 


65 


X2 


37.3 


-12.0 


24.9 


86.8 


108.9 


25 



TABLE II. Computed values of the self-energy corrections for 
modes Xi and X2, at zero temperature and zero frequency, 
and comparison of the renormalized frequency, Qz^, with the 
harmonic frequency Uu- All energies are in meV. The relative 
shift is seen to be quite large. 



coefficients of the form ^^J^,Vx(q) ^nd (q). 
These coefficients were obtained on a 4 x 4 x 4 q-mesh; 
Fourier interpolation was used to approximate them 
throughout the IBZ. 

The frequency dependent self-energy was computed in 
a range of interest, and the phonon spectral function 
was obtained^. In particular, by using the Lehmann 
representation^, it can be shown that 



(35) 



where 



but 3 to 4 times smaller than what the naive estimate 
based on dipsersionless parameters suggested. The cu- 
bic terms, which are real and negative at zero frequency, 
further reduce the estimate of the total self energy. The 
latter remains large enough however to strongly renor- 
malize with respect to cj^y, yielding a sizeable relative 
shift of 65 % ( 25 %) of the frequency of mode Xi (X2). 

The contributions to the linewidths coming from an- 
harmonic effects were also computed at and 300 K, 
for a fixed lattice geometry (a fixed value of a). Al- 
though a proper analysis of the temperature dependence 
should include the lattice expansion (which is beyond 
the scope of this work), we expect that the fixed cell 
calculations should still be indicative of the behavior of 
the linewidths at fixed pressure. The widths were found 
to be vanishingly small (0.3 meV) for mode Xi and 
0.6 meV (4.6 meV) for mode X2 at OK (300 K). Given 
that the electron-phonon contribution to the linewidth 
at X2 is about lOmeV, it should be possible to observe 
the temperature dependent contribution to this mode's 
linewidth, providing a possible experimental signature of 
the effects described here. 

In the approximation that anharmonicity does not af- 
fect the electron-phonon coupling (a reasonable assump- 
tion given that the g parameters are obtained from the 
deformation potential method, which is not affected by 
anharmonicity), the mode coupling is renormalized by 
anharmonicity and becomes 



\{anh) 



1 



TThN{0) l^<.(q,0)2' 



(37) 



A complete calculation of the renormalized value X^^^^) 
from equation (j4]) would require knowledge of the an- 
harmonic coefficients at points other than qx and is be- 
yond the scope of this work. However we can estimate 
an upper bound for the effect of the anharmonicity on 
the electron-phonon coupling. From Fig. [2] it is reason- 
able to assume that the electron-phonon coupling from 
to ~ 35 meV can be attributed to the region near Xi 
(partial value Ai 0^ 0.23) and that the electron-phonon 
coupling from ~ 85meV to ~ 135 me V can be attributed 
to the region near X2 (partial value A2 — 0.2); the rest of 
the coupHng is lumped together and assumed unaffected 
by anharmonicity (partial value A3 = A — Ai — A2 — 0.18). 
By estimating that the harmonic Eliashberg function is 
composed of two properly normalized S peaks at uji and 
UJ2 (the frequencies corresponding to Xi and X2) plus 
features unaffected by anharmonicity away from these 
frequencies, the renormalized coupling is given by 



(^)\.(g)\.A3^0.39, 



(38) 



(nnUq, io)) ' = (hio^iq)) ' + 2tkoMA^, '^)- (36) ^^^^^ ^"^^ ^2 are the renormalized frequencies of Xi 



The computed values of the "loop" and "bubble" self 
energy diagrams, at zero temperature and zero fre- 
quency, can be seen in Table [III along with the values 
of l]j,(qx,0). The "loop" contributions are quite large. 



and X2. A similar treatment yield ^cj^^ ^ 125 meV; 
using these renormalized parameters, the Allen-Dynes 
modification to the McMillan formula yields 1.5 K < Tc ^ 
5K for 0.14 > /i* > 0.1, suggesting that the anharmonic 
renormalization of the phonon spectrum, which acts to 
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FIG. 4. (Color online) Frozen phonon total energy calculations for the modes at X with frequency ^ 20 meV (Xi) and ~ 85 meV 
(X2), as a function of a displacement parameter 77. The values of the parameters huo and $4 are obtained from finite difference 
schemes. 



stiffen the modes which provides most of the contribution 
to A, leads to a strong reduction on Tc compared to the 
results obtained from the harmonic phonon spectrum. 



VI. CONCLUSION 

In this work, a full analysis of anharmonic effects to 
lowest order has been presented for modes in a region 
of the IBZ which provides most of the contribution to 
A, and this analysis has been compared to naive frozen 
phonon calculations. It has been shown that the modes 
Xi and X2 of AIH3 are strongly renormalized by anhar- 
monicity, although less so than a frozen phonon calcu- 
lation might have suggested. Indeed, it is found that 
the frequency of Xi is renormalized to 31.9 meV from 
19.4 meV (a 65% shift) and that the frequency of X2 is 
renormalized to 108.9 meV from 86.8 meV (a 25% shift). 
Furthermore, it is expected that anharmonicity induces 
a large, temperature dependent contribution to the X2 
mode linewidth, which could provide an experimental 
signature of the effect. A rough estimate suggests that 
renormalization could lead to a great reduction of the 
computed value of Tc, which could be as low as 2K, 
compared to the harmonic prediction of Tc— 20K. 

Thus, anharmonicity in AIH3 might play a role in 
explaining why the measured and computed supercon- 
ducting transition temperatures are in qualitative dis- 
agreement. There could be further phonon frequency 
renormalization due to the large electron-phonon cou- 
pling and, given that the electronic bandwidth near the 
Fermi energy is quite small, non-adiabatic effects in this 
system could also be substantial. 
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Appendix A: Symmetries of anharmonic coefficients 

The anharmonic coefficients can be defined as 

S^/7[{u}] 



$---:(Ri,...,RO 



du^l{Ki)...duZ{Kn) 



(Al) 

u=0 



which immediately implies that they are real and sym- 
metric under permutations of indices. Define the Fourier 
transformed coefficients as 



N 

{R} 

x$^-X"(Ri,...,R„).(A2) 
Translational symmetry, which can be expressed as 
$^--"(Ri,...,R„) = 

+ + (A3) 

for all lattice vectors R, implies that <l>^J ;^^(qi, ...,qn) 
vanishes unless qi + ... + is a reciprocal lattice vector. 
Furthermore, 

*".\-::"„"(qi, -.qn) = Kl'/Zi^i + ^i, -.qn + g„) 



(A4) 



for any set of reciprocal lattice vectors {G}. 



Appendix B: Extracting anharmonic parameters 
from dynamical matrices 
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The lowest order anharmonic correction to the self en- 
ergy, 

n = n(^) + n(^)+n(^), (bi) 
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involves anharmonic coefficients of third ("bubble" and 
"tadpole") and fourth ("loop") order. The most efficient 
and elegant way of obtaining these parameters is the use 
of the 2n + 1 theorem within the context of density- 
functional perturbation theory^. Briefly, this theorem 
guarantees that derivatives of the total energy up to order 
2n + 1 can be obtained from knowledge of the derivatives 
of the wave- functions to order n. In practice, however, 
only the ffist order derivatives of the wavefunctions {ie 
n = 1) are readily available, and finite-difference schemes 
are employed to obtain fourth order coefficients. 

An alternative way of obtaining the necessary coeffi- 
cients is through the frozen phonon method and finite dif- 
ferencing. What this method lacks in elegance, it makes 
up for by its simplicity and straightforward use, with- 
out the need for specialized software. The frozen phonon 
approach for this purpose is impractical for arbitrary q 
in the IBZ, as it implies computations with potentially 
very large super cells. However, the interesting point in 
this case is qx = 7r/a(0,0, 1), which lies on the side of 
the zone. 

Consider displacements of the ions from their equilib- 
rium positions of the form 



A5:(R;7?) 



.(qx), 



(B2) 



where 77 is a small real number and u is a specific mode 
of interest (Xi or X2). In Fourier space, this corresponds 
to 



to apply standard computational methods to extract the 
dynamical matrix. The displacements are periodic with 
respect to a lattice whose cells (henceforth supercell) con- 
tain two of the original cells stacked in the c direction. 
These subcells of the supercells will be labeled with ^ = 
and 5 = 1. Define 



A = a(0,0, 1); 



(B7) 



any lattice vector of the original lattice R can be ex- 
pressed as 



R = R + (5 A 



(B8) 



where R is some vector of the superlattice and S can be 
either or 1. The ionic positions with respect to this 
new basis are expressed as 

r(^,5)(R) = R + JA + + r/e'^'^x, + U(^,5)(R), 

(B9) 

where now the compound index (tz^S) identifies all the 
ions in the supercell with a label indicating ion and 
subcell 6. In reciprocal space. 



U(,,,) (R) = \/ ^ E (q) (BIO) 



A62(q;r?) = r?(5q,q^ViV<,(qx). 



(B3) 



The positions of the ions, now considered as simple num- 
bers and not operators, can be defined as 

r« (R) = R + b« + Ab« (R; r?) + (R) . (B4) 

The dynamical matrix computed about the non- 
equilibrium position is given by 



where q is a wave vector in the IBZ of the superlattice. 
The following relationships are thus immediate: 



0(^,5) (R) = u«,(R = R ,5A) 



U(«,<5)(q) = -y|-(^u«(q) + e' ''u«(q-FqjY)), 



(Bll) 



(B12) 



r,2 



C/[{u + Ab}] 



(B5) 

u=0 

(B6) 



from this last expression it is clear that the coefficients of 
interest for the computation of the "loop" diagram can 
be extracted from the second derivative with respect to 
77 of the out-of-equilibrium dynamical matrix. 

The situation is slightly more complicated, however. 
The ionic displacements introduced do not have the peri- 
odicity of the lattice; such periodicity is essential in order 



which implies 



-i5q- A 



Va<(q) 



d 



(B13) 



Above "q + qx" is meant to represent the wave vector 
inside the IBZ of the original lattice which is obtained 
from q + qx by an appropriate reciprocal lattice transla- 
tion. For what follows, it will be useful to remember that 
2qx is a reciprocal lattice vector of the original system, 
such that all anharmonic coefficients are periodic under 
q^q + 2qx. 

The dynamical matrix in the supercell representation 
is given by 



9 



r)aia2 (^\ 

^(/.,5)i,(/.,5)2 W 



d'U[{u}] 



gi(5i-(52)q-A 



u=0 



^".^".nq; V) + e^^'^+'^^-D^i:^ (q + q^; ^) + e^^^^O^^^"! (q; r?) + e^^^^O^^^^.^ (q + q^; r?) 



(B14) 



(B15) 



where 



and 



a<(-q)5^^"l(q + qx) 



a2[/[{u}] 



^ E *"^"l"3 (q- -q - qx, qx)x^3^(qx) + o(r?'), (bi6) 



- + ; ^) - an-(-q-qx)9n-(q) 



U=0 /^3,Q!3 



= ^ E *"!".'"3'(q + qx, -q,qx)<^(q^) + o{v^). (bit) 

u=0 /^3,a3 



From these last two terms, the anharmonic coefficients 
necessary for the computation of the "bubble" diagram 
can be extracted. 

Thus, from the supercell dynamical matrices D{r]) (the 
quantities actually computed through DFPT), it is a sim- 
ple matter of algebra to extract D{r]) and 0{r]). It is then 
useful to define partially projected anharmonic parame- 
ters, 

*"^:.Vx.x(q) ^ y E E <3^(q^)C(q^) 

X Kl^^:!:: (q, -q, qx , qx ) (B18) 
= |^^"^«|(q;^) + 0(r/), (B19) 

and 



These parameters can then straightforwardly be ob- 
tained by estimating the rj derivatives by finite differ- 
ence schemes. Furthermore, since these coefficients are 
periodic in reciprocal space (namely periodic under q 
q + G), they are susceptible to Fourier interpolation, in 
complete analogy with the usual procedures employed 
with dynamical matrices. Once these coefficients are ob- 
tained on a dense q mesh, it is possible to compute the 
"loop" and "bubble" diagrams. 



*"j".Vx(q) = ^EE<'-(q^) 

K3 as 

x*2™(q,-q-qx,qx) (B20) 
= -i-0:i:^i^;v)+Oiv). (B21) 
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